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Abstract — We propose a convex optimization procedure for 
identification of nonlinear systems tliat exhibit stable limit 
cycles. It extends the "robust identiflcation error" framework in 
which a convex upper bound on simulation error is optimized to 
fit rational polynomial models with a strong stability guarantee. 
In this work, we relax the stability constraint using the concepts 
of transverse dynamics and orbital stability, thus allowing 
systems vnth autonomous oscillations to be identified. The 
resulting optimization problem is convex, and an approximate 
simulation-error bound is proved without assuming that the 
true system is in the model class, or that the number of 
measurements goes to infinity. The method is illustrated by 
identifying a high-fidelity model from experimental recordings 
of a five rat hippocampal neuron in culture. 

I. Introduction 

Black-box identification of highly nonUnear systems poses 
many challenges: flexibility of representation, efficient opti- 
mization of parameters, model stability, and accurate long- 
term simulation fits, to name but a few [1], [2]. It is 
especially challenging when the system exhibits autonomous 
oscillations: such a system is intrinsically nonlinear and lives 
on the "edge of stability", since periodic solutions must have 
at least one critically-stable Lyapunov exponent [3]. 

Recently, a new framework has been introduced for 
identifying a broad class of nonlinear systems along with 
certificates of model stability and accuracy of long-term 
predictions [4]. However this method necessarily fails if the 
system has autonomous oscillations. In this paper we extend 
the method of [4] to remove this restriction. 

The main contribution of this paper is a method to identify 
highly nonlinear systems which: 

• searches over a very broad class of models, including 
those with limit cycles, 

• guarantees a (local) bound on deviation of open-loop 
model simulation from the data records, 

• is posed as a convex optimization problem, 

• is analysed without assuming that the true system is in 
the model class or that the number of measurements 
grows to infinity. 

A. Identification of Oscillating Systems 

In many scientific fields there is a need to capture oscilla- 
tory behaviour in the form of a compact mathematical model 
which can then be used for simulation, analysis, or control 
design. When the data comes from experimental recordings, 
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this is known as system identification. It is also becoming 
more frequent to perform model-order reduction via system 
identification methods from solutions of a very high dimen- 
sional simulation, e.g. computational fluid dynamics [5] or a 
detailed electronic circuit model (see, e.g., [6], [7]). 

In biology, systems that oscillate seem to be the rule 
rather than the exception: heartbeats, firefly synchronization, 
circadian rhythms, neuron spking, and many others [14], 
[15]. Nonlinear oscillator models have been used in speech 
analysis and synthesis, where stability of the identified model 
has been acknowledged as a major issue [16]. 

To the authors' knowledge, there is no generally appUcable 
methods of system identification - or model-order reduction 
- for oscillating systems. One family of approaches popular 
for aerospace model reduction is harmonic balance (describ- 
ing function) methods, in which the period of oscillation 
is assumed known and the model is reduced consider the 
problem in the fourier-series domain [8], [9], [5]. A similar 
approach has been taken to analyse phase-locked loops and 
oscillators, in which a local phase-offset system is of primary 
interest [7]. Neither of these approaches extend easily to 
situations in which the frequency of oscillation is input- 
dependent. Other papers assume a known decomposition into 
a stable Unear part and a static nonlinear map, and consider 
it a problem of closed-loop linear system identification 
[10]. Applications have included identification of combus- 
tion instabihties [11], [12]. A mixed empirical/physics-based 
approach has been used to produce low-order models of 
periodic vortex shedding in fluid flows [13]. 

B. Stability of Oscillations 

No linear system can produce an asymptotically stable 
limit cycle. Identifying nonlinear models from data is a dif- 
ficult problem, in particular because of the complex relation- 
ship between system parameters and long term behaviour of 
solutions. A recent approach, which this paper builds upon, 
works via convex optimization of a robust identification 
error which imposes an asymptotic stability constraint on 
the identified model [4]. 

However, if the system has a periodic solution, not driven 
by a periodic forcing term then this approach must fail: 
the stability constraint is too strong. To see this, suppose 
a system 

X = f{x) e R" 
has a non-trivial T-periodic solution x*{t), then x*{t -\- 



t),t G (0,T) is also a solution which will never converge 
to x*{t). 

The natural notion of stability for oscillating systems is 
orbital stability. A T-periodic solution x* is orbitally stable 
if nearby initial conditions coverge to the solution set in state 
space X = {x{t) : t E [0,T]} and not necessarily to the 
particular time solution x*{t). This is a weaker condition 
than standard (Lyapunov) asymptotic stability. 

Orbital stability can be studied via the introduction of so- 
called transverse coordinates, also referred to as the moving 
Poincare section [3], [17]. The basic idea is to contruct a new 
coordinate system at each point of the solution, decomposing 
the state into a scalar component tangential to the solution 
curve, and a component of dimension n — 1 transversal (often 
orthogonal) to the solution curve. 

It is known that periodic solution of a nonlinear differential 
equation is orbitally stable if and only the dynamics in 
the transverse coordinates are stable [3, Chap. VI]. This 
framework has previously been used to design stabilizing 
controllers and analyze regions of attraction for oscillating 
systems [18], [19], [20], [21], [22]. It has also been used to 
analyze the convergence of prediction-error methods when 
identifying a linear/static-nonlinearity feedback interconnec- 
tion that can oscillate [10]. In this chapter we extended the 
robust identification error method of [4] using a storage 
function in the transverse coordinates so as to robustly 
identify a broad class of nonlinear systems that may (or may 
not) admit autonomous oscillations. 

C. Paper outline 

The structure of the paper is as follows: in Section [ll] we 
set up the mathematical problems statement; in Section III 
we review the method proposed in [4] and explain why it is 



not suitable for oscillating systems; in Section IV we outline 
proposed approach and prove the main theoretical results; 
in Section |V] we give a convex (semidefinite) relaxation 
of the associated optimization problem; in Section VI we 



discuss practical matters of implementation and the utility 



of the model class; in Section VII we present experimental 



results fitting membrane potential dynamics of a spiking rat 
hippocampal neuron in culture; Section [VIII| has some brief 
conclusions; in two appendices we provide details of the 
experimental setup and some technical lemmas used in the 
proof of the main result. 

II. Problem Statement 

Given a data recorcQ of states, inputs, and outputs 
u(t), t e [0,r], the general problem is to 

construct a compact model in the form of a differential 
equation that, when simulated, faithfully reproduces the data. 
To pose the problem exactly we must specify both a model 
class to search over, and an optimization objective. 

'Here we assume that the data record consists of sufficiently smooth 
continuous-time signals on an interval, though in practice it will consist of 
a finite sequence of data points. 



A. Model Class 

The model class we will search over consists of 
continuous-time state-space models with state x E M", input 
u E W", output y E MP, and dynamics defined in the 
following implicit form: 



dt 



e{x) 



y 



f{x,u) 
g{x,u) 



(1) 
(2) 



where e : M" M",/ : M" x M™ : M" x 

M™ — W are smooth functions. The Jacobians with respect 
to X of e{x), f{x,u), and g{x,u) are denoted E{x) = 
■§^e{x), F{x,u) = ^f{x,u),G{x,u) = ^g{x,u). We 
will enforce the constraint that E{x) be nonsingular, so the 
above implicit model can equivalently be written in explicit 
form: 

X = E{x)~^f{x, u). 

Remark 1: To implement the methods described in this 
paper, e{x), f{x, u), and g{x, u) should come from a convex 
class of matrix functions for which we can efficiently check 
positivity. In practice, we use matrices of polynomials or 
trigonometric polynomials, so as to make use of sum-of- 
squares programming [23], [24]. 



B. Optimization Objective 

The general problem we consider is to minimize, over 
choice of e, f,g, the value of the simulation error: 



\y{t) ~ mi^dt 



where y{t) is the solution of ([TJ, (|2| with a;(0) = x{0). One 
may also wish to ensure that the dynamical system defined by 
([T]i, (|2| is well-posed and has some sort of stability property. 
Note that we do not assume that the system from which data 
is recorded is in the model class. 

Direct optimization of simulation error is not usually 
tractable: the relationship between system parameters and 
model simulation is highly nonlinear, and for black-box mod- 
els we typically don't have good initial parameter guesses. 
We make the problem tractable (a convex program) through 
a series of approximations and relaxations. 

III. Nonlinear System Identification via Robust 
Identification Error 

In this section we briefly present some results from [4] 
and explain why they cannot be directly applied to systems 
with autonomous oscillations. The basic idea is to search 
jointly for system equations as well as a storage function 
with output reproduction error as a supply rate. Standard 
dissipation inequality arguments [25] then provide a bound 
on long-term simulation error. 



A. Linearized Simulation Error 

Suppose we have a model of the form ([TJ, (|2]l and a 
data record {x(t),u{t),y{t)},t S [0,T]. We introduce the 
linearized simulation error as a local measure of the model's 
divergence from the data. 

First, we define the equation error signals associated with 
([T), ^ and the data: 

e^it) = E{i{t))i{t)~ f{i{t),ii{t)), (3) 
eyit) = y{t)-g{i{t),uit)). (4) 

Now, consider the following family of systems parametrized 

by 9 € [0, 1]: 

Eix)x = f{x,u) + fe, (5) 
y = 9{x,u)+ge. (6) 

Let {xe,ye) be the solution of the above system with fg = 
(1 — 9)€x and gg = {\ — 9)ey. That is, for 9 = 1 we have 
X0 = x,yg — y and for 6* = we have xg — x,yg = y. We 
can consider the following linearized simulation error about 
the recorded trajectory: 

f = lim -^ / \yg{t) - mi^dt 

as local approximation of the true simulation error £. 

B. Robust Identification Error 

Note that £ can alternately be represented as 

£ = / \G{x{t),u{t))A + ey{t)\^dt 
Jo 

where 

Ait) = \im -[xeit) - i{t)] (7) 
0~¥O a 

which obeys the dynamics 

j^E{i{t)W) = F{i{t)Mt))A{t) + e,(i). 

That is, A(i) is an estimate of the deviation of the model 
simulation x{t) from the recorded data trajectory x{t). 
It was shown in [4] that 

£< j £Q{t)dt (8) 
Jo 

for any g = g' > 0, wher^ 

£Q{t)= sup{2A'S'g(FA + e,) + |GA + ej,|2}. (9) 

The systems theory interpretation of (|9]l is that the first 
term in the supremum is the derivative of a positive-definite 
storage function with respect to linearized simulation error, 
and the second term is the output reproduction error. 

The bound ([HJ suggests searching over functions e, /, g 
and a matrix Q = Q' > so as to minimize the right- 
hand-side of (|8]l. This optimization is still non-convex, but a 

-Here, and frequently throughout the paper, we drop the arguments on 
E{x{t)),F(x{t),u(t)),G{x{t),u{t)),ex{t), and ty{t) for the sake of 
compactness of notation. It should be understood that these are always 
functions of time and the data. 



convex relaxation is given in [4] (we use a similar relaxation 
in Section [V] of the present paper). 

Each of the supremums over A in (|9]l are finite if and only 
if the matrices 

R = E'QF + F'QE + G'G 

for each data point is negative semidefinite. If this property 
holds for all x,u, then it has been proven that the system 
is globally incrementally output stable. Roughly speaking: 
A'E'QEA is a Lyapunov function for the difference be- 
tween two solutions A « xi-X2, and A'{E'QF+F'QE)A 
is its derivative. A formal proof of stability is given in [4]. 

For the purposes of the present paper, it is sufficient to 
note that enforcing global incremental stability is too strong 
to allow identification of systems exhibiting autonomous 
oscillations, since such systems cannot satisfy this property. 
The main purpose of this paper is to overcome this limitation 
via a reformulation of the RIE in the transverse dynamics. 

IV. Transverse Robust Identification Error 

In oscillating systems, perturbations in phase cannot be 
stable and will therefore accumulate over time. The natural 
form of stability is orbital stability, which can be defined 
as stability to a solution set in state space, rather than a 
particular time solution. A standard framework for anlaysis 
of orbital stability is via transverse coordinates. 

Correspondingly, if both a true system and an identified 
model admit autonomous oscillations, then it is not possible 
to ensure that phase deviations between them converge in 
time. Therefore we introduce the concept of orbital simula- 
tion error: 

r \y{t) - yirm'dt 
Jo 

where 

T(t) — arg min \x(t) — x(t)\. 
re[o,T] 

We also define the approximation linearized orbital sim- 
ulation error. 

£:^:=liml^ \yeit)-y{Tg{t))fdt 

where 

Tg{t) = arg min \xg{t) - x{t)\ 
Te[o.T] 

Observe that To{t) — t and Te{t) is well-defined and smooth 
function of 9 when 9 is sufficiently close to zero. We will 
never explicitly compute Tg. 
This can be rewritten as 

£^ := r\GA{t) + ey\' 
Jo 

where 

A{t) = \im -[xg{t) - iiTg{t))]. (10) 

Remark 2: The difference between A in(|7]) and A in ( [T0| 
is as follows: A{t) is an estimate of the deviation of the 
model simulation x{t) from the recorded data trajectory x{t) 



at the same time t; whereas A{t) is an estimate of the 
deviation of the model simulation x{t) from the closest point 
on the data trajectory x{t) where r = argminsg[o t] \x{t) — 
xis)\. 

A. Simulation Error Bound 

We are now ready to state the main theoretical result of 
the paper, which gives an upper bound for the linearized 
simulation error in terms of a point-wise condition on the 
data and model parameters. 

Define the following projection operators: 

_ x{t)x{t)' :=I-n{t), (11) 

i.e. Tr{t) projects on to the one-dimeonsional subspace par- 
allel to x{t) and Ii{t) projects on to the {n — 1) -dimensional 
subspace orthogonal to this. 
Let 

£^{t) = sup {2An'^'£;'Q((i^ + £;ri)n'^A + e^) 



\GWA + ey\'} 



(12) 



where 11'" (i) € is a matrix with orthonormal rows 

spanning the subspace orthogonal to x and Q a symmetric 
positive-definite n x n matrix. It is a "reduced" form of the 
rank {n — 1) matrix Il{t) containing only independent rows. 

Theorem 1: Consider measurement x{t),y{t) and simula- 
tion x{t),y{t) with the same initial condition i(0) = x{0) 
and the same input u{t). Then the following relation holds: 



lo 



(13) 



Proof: The inequality is shown via a dissipation inequality 
for the following storage function: 

y(A, t) = \Eii{t))n{t)A\l + \n{t)A\\ (14) 

Consider in particular, for each time t, the value V{A{t), t). 
Note that, due to the particular structure of the storage 
function, since A is infinitesimal, 

V{A{t),t) = min ViAM.t) 

■re[0,T] 

Using Lemma [2] (see appendix) 

jV{A(t),t) + \GA{t) + ey\'<£^{t) (15) 

Since 5(0) = a;(0) we have V^(A(0),0) = 0, so integrat- 
ing ([TSj gives 

y(A(T),T)+ / \GA{t)+ey\^dt< [ £{t)dt. 
Jo Jo 

and by definition V{A{T),T) > 0, so the above inequality 
implies ( pj] ). This completes the proof of the theorem. 

□ 



V. A Convex Upper Bound 
Theorem [T] suggests minimizing 

over choices of e,f,g and Q as an effective procedure 
for system identification. However, this is still a nonconvex 
optimization. In this section we propose a convex upper 
bound for which one can efficiently find the global minimum 
via semidefinite programming. 

The basic idea is to decompose the each non-convex term 
into the sum of a convex and a concave part, and upper- 
bound the concave part with a linear relaxation. 

Theorem 2: Define 

A+ = Eii)iI + U)WA + FiS:,u)U'^A + e, 
A- = E{i) {I - t{)WA -F{i, u)WA - 



GIT A + Ey. 



Then E^it) < £^{t) where 



£q = sup 



- (n'-A)'A- + |A, 



which is convex in e, /, g, and Q^^ > 0. ■ 

Proof: A similar statement was proved in [4, Section V]. 
The upper bound is based on the expansion 

\a - Q"^A|| = AQ^^A - 2A'a + a'Qa 

which, when Q > 0, clearly implies 



-a'Qa< A'Q-^A-2A'a. 



(16) 



Notice that the right-hand side of ([16]) is convex in a and 
Q^^ whereas the left-hand-side is concave. 

Note that we also have the following expansion: 

AiEIl^AyQ[{F + Etl)WA + e,] = |A+|2,-|A-||. (17) 

The first term on the right-hand-side of ([TtJ is convex in e, /, 
and Q^^ and the second term is concave. Setting a = A~ 
and applying ([T6| to ^TT) gives the statement of the theorem. 

□ 

Summarizing the results of this and the previous sections, 
we have the relations 

(£^«£^< r £^{t)dt< r £^{t)dt, 

Jo Jo 
with the leftmost term being the true orbital simulation 
error, and the rightmost term being convex in the system 
equations and Q. Therefore we propose to perform system 
identification via the optimization 

i-T 

£Q(t)dt — > min 

over choices of e, /, g, and Q^^ > subject to the constraint 
E{xy + E{x) > I for all x. 



In practice, we will have a record of the true system at a 
finite number of times ti,i — 1,2, ...,N, and as a surrogate 
for the above we minimize the finite sum of the TRIE terms: 

N 

'^S^iU) min. 

1=1 

VI. Implementation and Discussion 

We now discuss practical considerations for data prepara- 
tion and minimization of the upper bounds using semidefinite 
programming. 

A. Extracting States from Input/Output Data 

The RIE formulation assumes access to approximate state 
observations, x{t). In most cases of interest, the full state of 
the system is not directly measurable and extraction of a state 
vector is a challenging problem in its own right. In practice, 
our solutions have been motivated by the assumption that 
future output can be approximated as a function of recent 
input-output history and future input. A common method 
for choosing a set of time-delays is to optimize mutual 
information [26]. To summarize recent history, we have had 
success applying linear filter banks, as is common in linear 
identification (e.g. Laguerre filters [27]). 

Projection-based methods such as subspace identification 
[28] and proper-orthogonal-decomposition [29] are common 
methods for linear system identification and model reduction. 
However, even in fairly benign cases one expects the input- 
output histories to live near a nonlinear submanifold of the 
space of possible histories. Connections between nonlinear 
dimensionality reduction and system identification have been 
explored in some recent papers, e.g. [30] and [31]. 

For experimental recordings, derivatives can be estimated 
via differentiation filters or noncausal smoothing before 
numerical differentiation. Approximating additional states 
through filter banks allows the rates of these variables to 
be calculated analytically. 

B. Semidefinite/Sum-of-Squares Programming Formulation 

For each data point, the upper bound on the local TRIE is 
the supremum of a concave quadratic form in A. So long as 
e, / and g are chosen to be linear in the decision variables, 
this upper bound can be minimized by introducing an LMI 
for each data-point; we introduce a slack variable Si for each 
data-point: 

S^>£Q{ti), 

which can be transformed via Schur complement into an 
LMI constraint. Then we optimize for Si min. We 
parametrize e, /, and g as polynomials, so that a sum-of- 
squares relaxation [23] is used to enforce the well-posedness 
consti-aint E{x) + E{x)' > I \/x e M". 

Note that the robust identification criterion typically has 
a relatively small number of decision variables, dependent 
on the order and degrees of the system model. However, to 
transform the problem into a standard semidefinite pogram, 
a very large number of slack variables, equality constraints, 
and LMI constraints are introduced, growing with the number 



of data points. The main reason for doing this is to make 
use of existing and well- tested semidefinite solvers, however 
it is possible that other methods - such as cutting-plane 
algorithms - will be more computationally efficient (see, e.g., 
[32]). This will be a topic of future work. 

C. Scientific Utility of the Model Class 

A disadvantage of our model class is that parameters have 
no physical meaning. On the other hand, direct identifica- 
tion of physics-based models can have many local minima 
[33] or be very ill-conditioned (i.e. the models are weakly 
identifiable) [34]. In such cases, it is an acknowledged issue 
that scientific inference from identified physical parameters 
could be difficult or misleading. 

Whilst the parameters of our model class have little 
meaning on their own, this disadvantage is offset somewhat 
by the fact that the models themselves - being polynomial 
in form - are well-suited to analysis via sum-of-squares 
programming, an advantage not shared by other broad classes 
of nonlinear system models. For example, one can estimate 
regions of attraction of limit cycles [20], perform worst-case 
or stochastic safety verification [35], and model validation 
[36]. These methods have recently seen growing application 
in systems biology [37] [38]. 

VII. Experimental Results on Live Neurons 

We now demonstrate the method by identifying the 
membrane-potential dynamics of a live, in-vitro hippocampal 
neuron. A micropipette is used to establish an interface with 
the cell such that current can be injected into the soma, 
and the membrane potential response can be recorded. A 
microscopic photograph of the neuron and patch is shown in 
Figure [T] The preparation of the culture is described in the 
appendix. 

The membrane dynamics of a single neuron are highly 
complex: at low currents the system responds like a low- 
order passive linear system. However, when certain ion 
channels are activated rapid spiking can occur After spiking 
there is a refractory period in which sensitivity is reduced, 
and with sufficient current input spiking can repeat at an 
input-dependent frequency. 

There is a spectrum of models of neuron dynamics, 
ranging from simple "integrate and fire" models to highly 
complex biophysical models of ion channels and conduc- 
tances [39]. Threshold based models generally have a very 
small number of parameters, but do not provide high fi- 
delity reproduction of the membrane potential dynamics. 
By contrast biophysical models can be very accurate, but 
are highly nonlinear and are very difficult to identify [40] 
- they can have many locally optimal fits in disconnected 
regions of parameter space [33]. In this section we use the 
proposed method to identify a black-box nonlinear model 
with comparatively few states (three) which reproduces the 
experimentally observed spiking and subthreshold behavior 
with very high fidelity. 

Three increasing step currents are applied to the neuron 
resulting in increasing firing rate and a characteristic change 
in the spike ampHtude and shape. 



100- ■ - 

200 ; # 

300 - 

400 - - 

500 - . tBIIP 

600 • 

700 - ^ * ^ - 

800 - 

900 - .... - 

1000 - 

200 400 600 BOO 1000 1200 

Fig. 1. The neuron in culture and the glass micropipette electrode used to 
interface to it. Imaging: phase contrast image at 20 X magnification on an 
inverted Olympus IX-71 microscope. Scale: 100 pixels (marked on axes) = 
43 microns. 

As discussed in Section |VI] we must find a good proxy for 
the internal state of the system. Here we used two Laguerre 
filters with identical pole locations to summarize the recent 
history voltage history. 

Figure [2] presents a comparison of fit performance using 
three methods. The first is equation error minimization, i.e. 
simply optimizing 

^ \E(x{U))5:{U) - j{i{U),u{U))\^ ^ min 

i 

subject to the well-posedness constraint E{x) + E{x)' > I 
but without constraints on stability or long-term simulation 
error (this is similar in principle to NARX and prediction 
error methods). The second method is the comparison is 
the original RIE minimization from [4], and the third is the 
proposed Transverse RIE method. 

We see that while equation error minimization (top) leads 
to initially good performance, the model goes unstable quite 
quickly. Fitting with the RIE (middle) leads to the anticipated 
overly stable model dynamics (see Section|lIl]). The final plot 
presents the Transverse RIE identification, which matches the 
experimentally observed spike patterns very well. 

We have also had success identifying behavior which 
covers both the subthreshold and spiking regime of a neuron. 
The applied stimulus was a variety of multisine signals. 
Figure [3] presents validation of a Transverse RIE fit on held- 
out data. The lower plot is the multisine input in pico- 
Amperes. The upper plot presents the original data and fit. 
Both the subthreshold regime and spikes are generally well 
reproduced. 

VIII. Conclusions 

This paper has introduced a new technique for identifica- 
tion of nonlinear systems which may produce autonomous 
oscillations, i.e. system oscillations which are produced 
internally by the dynamics rather than as a response to a 
periodic input. 




Fig. 2. A neuron was subjected to increasing step-cuiTcnts, and spiked 
with increasing frequency. Long-term simulation of an equation error fit 
(top) is unstable. RIE (middle) minimization provides an overly stable fit. 
The proposed TRIE method (bottom) reproduces the spikes accurately. 

A convex optimization procedure is developed which 
minimizes an upper bound on a local measure of simulation 
error - the long-term divergence between a model simulation 
and the recorded data. 

The proposed method worked well on the challenging 
problem of accurately modeling the membrane dynamics of 
a live neuron from experimental current-voltage recordings. 
The input-dependent absence or presence, and frequency of 
repitition, of spiking events was well captured in the model. 

Future work will include investigating dedicated solvers 
and methods for extracting states, as well as testing the 
proposed method on a wider range of systems, such as 
electronic oscillators and vortex-shedding fluid flows. 

Appendix 
A. Live Neuron Experimental Procedure 

Primary rat hippocampal cultures were prepared from PI 
rat pups, in accordance with the MIT Committee on Animal 
Care policies for the humane treatment of animals. Dissec- 
tion and dissociation of rat hippocampi were performed in a 
similar fashion to [41]. Dissociated neurons were plated at a 
density of 200K cells/mL on 12 mm round glass coverslips 
coated with 0.5 mg/mL rat tail collagen I (BD Biosciences) 
and 4 /ig/mL poly-D-lysine (Sigma) in 24-well plates. After 
2 days, 20 /iM Ara-C (Sigma) was added to prevent further 
growth of glia. 

Cultures were used for patch clamp recording after 10 days 
in vitro. Patch recording solutions were previously described 
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Fig. 3. Response of real neuron and TRIE model simulation, showing both 
stable sub-threshold behavior and spiking. 



in [42]. Glass pipette electrode resistance ranged from 2- 
4 Mrt. Recordings were established by forming a Gil seal 
between the tip of the pipette and the neuron membrane. 
Perforation of the neuron membrane by amphotericin-B (300 
/ig/mL) typically occurred within 5 minutes, with resulting 
access resistance in the range of 10-20 Mfi. Recordings with 
leak currents smaller than -100 pA were selected for analysis. 
Leak current was measured as the current required to voltage 
clamp the neuron at -70 mV. Synaptic activity was blocked 
with the addition of 10 ^iM CNQX, 100 APV, and 10 
bicuculline to the bath saline. Holding current was applied 
as necessary to compensate for leak current. 

B. Two Technical Lemmas 

Lemma 1: Given the definitions of A in (|7]), A in ( [TO] i, 

and n(i) in ([TT), thenfor all t e [0,T]: K{t) = n(i)A(i) 

Proof: By definition, 



where 



Tg = arg min lifr) — xg\ 

re[0,T] 



First note that for any 9 sufficiently small, 

X0{t) ~ xir) xo{t) — x{t) 



n(r)- 



e 



since x{t) is the closest point to the solution set to 
xo{t), and so xg{t) — x{t) is in the subspace orthogonal to 
x{t) and is therefore invariant under the projection n(r). 
And since n(T) — > n(t) as 6* ^ we have 

Ht) = W) liniJ(xeW-i(Te)). 

Now, define xo{t) = x{tq) — x{t). Since Tg — t — 0(6), 
it follows that xg{t)/d is a rescaling of x for 6 small, i.e. 

r xe{t) • 
lim — - — = axlt) 

9^0 9 



for some a e ]R. Now, by definition Il{t)x{t) = so 

A{t) = n{t) lim hxg{t) ~ i(t)) = n(i)A(i) 

□ 

Lemma 2: Given the storage function ( [T4] i, the following 
inequality holds for all t e [0,T]: 

jVCm^t) + |GA(t) + < E^(t) (18) 

where Eqit) is given in ( [T2) l. 
Proof: Let 

V = E{i{t))Tl{t)\ 5 = 7r(t)A, 

so V^(A,t) = v'Qv + 5'5 and 

V = 2v'Qv + 255 

but from Lemma [l] for the particular case A{t) ~ A(t) we 
have n(i)A(i) = A(i) and 7r(i)A(<) = 0, so 

V = v'Qi. 

We now derive an expression for v. 

Let ti ~ t and decompose A(t) into two components 
Aj^(t) + 5{t), transversal and tangential, with respect to 

x{h)- 

Note that this decomposition is based on the transversal and 
tangential decomposition at a. fixed time ti, not based on a 
rotating coordinate system. By the chain rule, 

[E{x{t))mm)] 



dt 
d_ 
dt 
d_ 
dt 



[E{x{t))\l{ti)A{t) + ^(x(ti))n(i)A(ti)] 

[i?(x(t))n(ti)(A^(<) + 5(<))] 



t=ti 



+S(a;(ti))n(ti)A(ti) 

but by definition IV{ti)6{t) = and n(ii)A_L(<) = A_L(i) 
for all t, so 

d 



dt 



[E{x{t))A^{t)] 



E{x{ti))n{ti)A{ti) 



= F{xiti),uiti))A^{ti) + + ^(.T(ii))ri(ii)A(ti) 

But Aj_{ti) = A{ti) = n(ii)A(ti), so for this particular 
A, 

V{A{t), t) = 2AIl'E'Q{{F + Etl)IlA + e^) (19) 

Now, since A{ti) = A{ti) and V{A{t),t) < V{A(t),t) for 
all t we have 

T>(A(ti),ii) < ■l/(A(ti),ti) (20) 
so given ([T9| and ( |20l ), 

y(A(ii),ti) + |GA(t)+e,p<f^(t) 
which proves the lemma. □ 
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